Supplementary Content

Functional diversity … Predictions of a model of an oxic-anoxic microbial ecosystem

1 Introduction

This document contains the Supporting Content for the report ?????????????.

2 Creating among strain diversity

Here we provide additional description of the method used to create strain variation within functional groups. The following should be read only after the relevant sections of the main text.

Take for example the cyanobacteria functional group, a reference maximum growth rate parameter \(G_{max,CB}\) and a reference tolerance parameter \(h_{SR,CB}\). Recall that amount of among strain variation is controlled by the two meta-parameters \(Div_{G_{max,CB}}\) and \(Div_{h_{SR,CB}}\). To create a tradeoff, the two meta-parameters were always of different sign. The among strain variation was calculated such that the growth rate of strain \(i = 1\) was a factor \(1/({2Div_{G_{max,CB}}})\) of the reference growth rate, and the growth rate of the \(i = n\) strain was a factor \(2Div_{G_{max,CB}}\) of the reference growth rate. The growth rate of the \(i = {2,...n-1}\) other strains was the reference growth rate multiplied by a factor that was equally distributed between the factor of the \(i = 1\) and \(i = n\) strain. This implementation of the variation and trade-off has the assumption of a linear trade-off of log-log transformed among-strain variation. Hence, for example, with the reference growth rate \(G_{max,CB} = 0.05\), among strain variation of \(Div_{G_{max,CB}} = 1\), and nine strains (\(n = 9\)), the growth rates of the nine strains are 0.025, 0.03, 0.035, 0.042, 0.05, 0.059, 0.071, 0.084, 0.1. Put another way, the growth rate parameter of the \(i\)-th of the \(n\) strains in a functional group was calculated as \(G_{max,i} = D_i * G_{max}\) where \(D_i\) is the \(i\)-th element of the set \(D = {2f(x)Div_{G_{max}} |x ∈ N, x ≤ n}\), where \(f(x) = (x − (n + 1)/2)/((n − 1)/2)\).

This procedure ensures that the range of trait values is independent of the number of strains. When there were only two strains, the two strains took the minimum and maximum traits values, as defined by the \(Div\) variable for that trait. When there were three strains, one had the minimum trait value, one the maximum, and one the mean. With nine strains, one had the minimum trait value, one the maximum, one the mean, and three were equally spaced (on log scale) between the minimum and the mean, and the remaining three between the maximum and the mean.

We show three examples of the created variation among strains within functional groups, one for nine strains (Figure 2.1), one for three (Figure 2.2), and one for two (Figure 2.3). In each, three graphs shows the diversity among strains in the maximum growth rate trait and the tolerance trait, plotted against each other and therefore also showing the tradeoff between the two traits.

Notice that in all three cases, the range of trait variation represented by the strains is the same, even though the number of strains differs (a consequence is that the lower row of panels in Figures 2.1, 2.2, and 2.3 are identical). This approach allowed us to independently manipulate trait diversity and number of strains.

All three examples have the following amounts of trait variation: \(Div_{G_{max,CB}}\) = 0.032, \(Div_{h_{SR,CB}}\) = -0.16, \(Div_{G_{max,SB}}\) = 0.095, \(Div_{h_{O,SB}}\) = -1.938, \(Div_{G_{max,PB}}\) = 0.095, \(Div_{h_{O,PB}}\) = -1.938. These are also the maximum amounts of trait variation explored in the simulations.

Figure 2.1: Trait variation with nine strains in each of the functional groups. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains.

Figure 2.2: Trait variation with three strains in each of the functional groups. Note that the three displayed strains are the 1st, 5th, and 9th strains from the nine strain example shown in Figure 2.1. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains.

Figure 2.3: Trait variation with two strains in each of the functional groups. Note that the two displayed strains are the 1st and 9th strains from the nine strain example shown in Figure 2.1. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains.

We also show, in Figure 2.4, the relationship between realised growth rate and the concentration of the respective inhibiting substance, e.g. sulfide for the cyanobacteria. In these graphs only the least and most tolerant strains are shown. The lines cross due to the tradeoff. At high concentrations of the inhibiting substance the more tolerant strain has higher realised growth rate, while at lower concentrations it is the less tolerant strain that has the higher realised growth rate.

Figure 2.4: Response of realised growth rate to variation in concentration of the inhibiting substrate for two strains, namely the most and the least tolerant strain. As always colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains.

3 Method for finding stable states

We now discuss two approaches for finding how the stable state of the system varies along the gradient of oxygen diffusivity when there exists the possibility of multiple stable states. In the main article we describe and report results of the method that we term the Temporal method.

Here we also describe and give results of what we term the Replication method. In this, we made a separate simulations for each of many values of oxygen diffusivity, ran the simulations for long enough to ensure transients had passed, and recorded the state. We ran multiple simulations per value of oxygen diffusivity, and started them with different initial conditions (either low or high total cyanobacterial abundance), so as to check for presence of alternate stable states. This method used by Bush et al 2017. We term it the replication method since it is achieved by having many independent replicates of the ecosystem, each with a different oxygen diffusivity, and also a separate replicate for each of the different chosen initial conditions.

As is shown in the figures below, the temporal method results in a larger region of bistability than the replication method, but only during increasing oxygen diffusivities. The difference in results occurs because, as mentioned above, the initial conditions differ between the two approaches. For example, in the temporal method, the oxygen concentration in the anoxic state is lower than the initial oxygen concentration in the replication approach, and Owen guesses that the sulfide concentration is higher check that last statement.

3.1 Replication method, no diversity

In the following graph, the solid (down) line refers to the stable state when cyanobacteria begin at high abundance (favouring the oxic state). The dashed up line refers to when their initial abundnace is low (favouring the anoxic state).Stable states found with the replication method. Please note: i) that near vertical lines indicate the position of a tipping point and hence the system cannot exist on point of that part of the line; ii) the separatrix (the unstable equilibrium) is not shown.

3.2 Temporal method, no diversity

Same as for the previous graph.

3.3 Both methods overlaid, no diversity

Stable states of the two methods overlaid, the temporal method with a dashed line, and the replication method with a solid line.

3.4 Effect of diversity, replication method

The replication method shows the same effects of functional diversity as those revealed by the temporal method:

• Introducing functional diversity in all three functional groups increases the region of bistability relative to when there is no diversity within all functional groups).
• Introducing functional diversity in only the cyanobacteria function group delays (with respect to the oxygen diffusivity it occurs at) the switch to the anoxic state as oxygen diffusivity is reduced, but has no effect on the oxygen diffusivity at which the anoxic to oxic transition occurs .
• Introducing function diversity in only the sulfur bacteria delays (with respect to the oxygen diffusivity it occurs at) the switch from anoxic to oxic, and advances the switch from oxic to anoxic.

These results are shown in the figure below, which shows results of the replication method.

The figure above shows the effect of response trait diversity (y-axis) on the positions of the tipping points (show by dark lines), which determine the extent of the region of bistability in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations.

3.5 Further thoughts

Here we give some further thoughts regarding the comparison of the temporal and replication method for stable state finding.

A particularly notable difference between the two approaches is in the initial condition. In the temporal approach, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the approach of Bush et al (2017) the initial conditions are fixed and arbitrary (though this approach is still useful for illustrating the bistability of the system). This difference in initial conditions increases the region of bistability in the temporal approach because as a tipping point approaches, the functional group that will dominate after the regime shift starts from much lower abundances than in the approach by Bush et al (2017) and thus requires stronger amelioration of the environment to take over. Nevertheless, our main conclusions about the effects of functional trait diversity are the same using either approach.

The temporal approach is particularly useful when investigating the effects of strain richness. In preliminary simulations using the Bush et al (2017) approach, we found the counterintuitive result that higher strain richness led to lower resilience. We then recalled that we chose to divide the initial total abundance of the cyanobacterial (CB) functional group equally among the strains (termed a substitutive design in biodiversity manipulation experiments). This meant that the most tolerant CB strain in a three-strain simulation had higher initial abundance than the same most tolerant CB strain in a nine-strain simulation. The greater abundance of the most CB tolerant strain made it possible for a three-strain system to develop to be oxic, whereas the nine-strain system developed to be anoxic, despite identical range of trait variation. We could have chosen to do otherwise, for example start with the least tolerant strains of cyanobacteria as more abundant when the system started with high cyanobacterial abundance, and the most tolerant strain as more abundant when the system starts with low cyanobacterial abundance. Using the temporal approach removed the need to specify initial abundances (except the very first one), and therefore results were not so determined by the assumptions used to determine the initial conditions.

4 Effect of simulation duration

Review what the wait time is.

Wait times of 1e6 hours produce the following results (same as those presented in the main report).

all_stab <- readRDS(here("experiments/0_ss_finding/temporal_method/processed_data/stab_data_temporal_method.RDS"))

wait_time <- 1e6
num_strains <- 9
plot_here <- all_stab %>%
  filter(waittime == wait_time)
p1 <- fig_div_vs_o2diff_1strain_7row(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(wait_time,
                                                           " wait time,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ") 
p1 + p2

wait_time <- 1e5
num_strains <- 9
plot_here <- all_stab %>%
  filter(waittime == wait_time)
p1 <- fig_div_vs_o2diff_1strain_7row(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(wait_time,
                                                           " wait time,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ") 
p1 + p2

wait_time <- 1e4
num_strains <- 9
plot_here <- all_stab %>%
  filter(waittime == wait_time)
p1 <- fig_div_vs_o2diff_1strain_7row(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(wait_time,
                                                           " wait time,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ") 
p1 + p2

wait_time <- 1e3
num_strains <- 9
plot_here <- all_stab %>%
  filter(waittime == wait_time)
p1 <- fig_div_vs_o2diff_1strain_7row(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(wait_time,
                                                           " wait time,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ") 
p1 + p2

wait_time <- 1e2
num_strains <- 9
plot_here <- all_stab %>%
  filter(waittime == wait_time)
p1 <- fig_div_vs_o2diff_1strain_7row(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(wait_time,
                                                           " wait time,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ") 
p1 + p2

Wait times of 1e5 produce the following results, which qualitatively and often quantitatively very similar to those for wait times of 1e6. There is a notable difference when there is variation in only SB, or CB-SB variation treatments with two or three strains:

Here, moderate increases in diversity have the same effect as with wait times of 1e6. However, larger increases in diversity suddenly cause a large reduction in resilience of the anoxic-oxic transition. Here is an example of the stable states when diversity is moderate:

wait_time <- "1e+04"
ss_3s <- readRDS(here(paste0("experiments/0_ss_finding/temporal_method/data/ss_data_3strains_waittime", wait_time, "_event_definition_2.RDS")))


#max <- max(ss_3s$PB_var_gmax_s)
#f <- 0.625
#fv <- f*max

only_SB_variable <- ss_3s %>%
  filter(
    abs(CB_var_gmax_s) < 0.0001,
    abs(PB_var_gmax_s) < 0.0001
  )
only_SB_variable <- only_SB_variable[order(only_SB_variable$SB_var_gmax_s), ]

only_SB_variable$SB_var_gmax_s

##  [1] 0.000000000 0.004986150 0.009972299 0.014958449 0.019944599 0.024930748
##  [7] 0.029916898 0.034903048 0.039889197 0.044875347 0.049861497 0.054847647
## [13] 0.059833796 0.064819946 0.069806096 0.074792245 0.079778395 0.084764545
## [19] 0.089750694 0.094736844

fig_state_vs_o2diff(only_SB_variable[13,])

Increases in oxygen diffusivity cause strain replacement in the SB to the more tolerant strains, which then delay the shift. When there is just a little more diversity in the SB, this happens:

fig_state_vs_o2diff(only_SB_variable[14,])

We believe that this phenomenon is caused by an effect that we explain with reference to an anology of using stepping stones to cross a river. Increasing diversity in our simulations is akin to increasing the width of the river, but while keeping the number of stepping stones (strains) constant. So long as the river is not too wide, we can still cross it, but once the width of the river is such that the stepping stones are too far apart, we can’t cross the river. And this happens suddenly, the width increases just a tiny bit more, and we suddenly can’t cross the river when we could before.

In the simulated microbial ecosystem, the strains are the stepping stones, and it is much hard (takes longer) for strains to replace each other when there is greater distance (in trait space) between them. When they can no longer replace each other in the time available, the system behaves as if only the least tolerant, highest growth rate strain (lighted shade of colour) is present. In the case illustrated above, this leads to much lower resilience of the anoxic system state.

This phenomena is dependent on having static trait values for the strains. If strain traits were dynamic (e.g. due to evolution), then we expect this result to not occur, and rather for the effects of diversity with two or three strains to be the same as those with six or nine.

4.0.1 Example of transient dynamics at 1’000’000 hours

The results below show the stable states as a function of oxygen diffusivity for a system with intermediate level of diversity in all functional groups, with either 1’000’000 hours simulation time, or 2’000’000 hours. This shows that except very close to the values of oxygen diffusitivity at which a state transition occurs, 1’000’000 hours duration of the simulations is sufficient to reach stable state, or very close to it. Very close to the transitions, even at 2’000’000 hours, the dynamics are still transient, with high tolerance strains at high density, even though they will eventually be outcompeted by low tolerance (higher growth rate) strains.

Below is Figure 3 from the main manuscript, in which the duration at each oxygen diffusivity value is 1e6:

For comparison, here are the results with twice the duration at each oxygen diffusivity:

5 Finding the separatrix

In the original publication (Bush et al 2017) Figure 3a shows a separatrix in the cyanobacteria abundance - oxygen diffusivity dimensions of the system. If the system has initial cyanobacterial abundance greater than the value at the separatrix then the cyanobacterial abundance increases, and the system goes to the oxic stable state. Conversely, if the system has initial cyanobacterial abundance greater than the value at the separatrix then the cyanobacterial abundance decreases, and the system goes to the anoxic stable state.

It is critical to note, however, that the position of the separatrix depends on both the parameters of the system, and the initial conditions of variables other than cyanobacterial abundance. If the system were started with higher initial concentration, then the position of the separatrix would change (as would the values of oxygen diffusivity at which the system tipped from anoxic to oxic, and from oxic to anoxic).

This issue of dependence of the position of the separatrix on the value of initial condition of multiple state variables has particular significance when there are multiple strains per functional group. Then one cannot just vary initial total cyanobacterial abundance (as was done to create the separatrix on Figure 3a in the original publication) but one must also define how the initial abundance of each of the strains varies. One might choose to make each strain have equal initial abundance, and to vary the total (as we did in our implementation of the replication method in the main article).

We could, however, and arguably should, give the most tolerant strain of cyanobacteria higher initial relative abundance when total initial cyanobacterial abundance was low, and the least tolerant strain higher initial relative abundance when total cyanobacterial abundance was high. And we would need to specify the total and relative abundance of all cyanobacteria (and sulfur bactria) strains. Thus, with a system as complex as the 9 strain system, which has 31 state variables the separatrix is very high dimensional, and evaluating (and plotting) it in one dimension requires many assumption. We chose to not attempt to find a meaningful set of assumptions in this study, though not because we think it would be unprofitable to do so, but rather because we did not want to risk distracting attention from findings we consider to be significant enough in the own right.